setup;
global h

%sim_u = 7*sim_x^7*sim_y^7+2*sim_x^3*sim_y^3+3*sim_x^9*sim_y^2+8*sim_x^6*sim_y^2;
%sim_f = diff(sim_u,sim_x,2) + diff(sim_u,sim_y,2);

x = 2 : h : 14;
y = 3 : h : 16;
[xx,yy] = meshgrid(x,y);

m = length(x);
n = length(y);

fu = @(x,y) 7.*x.^7.*y.^7+2.*x.^3.*y.^3+3.*x.^9.*y.^2+8.*x.^6.*y.^2;
ff = @(x,y) 294.*x.^5.*y.^7+294.*x.^7.*y.^5+12.*x.*y.^3+12.*x.^3.*y+216.*x.^7.*y.^2+240.*x.^4.*y.^2+6.*x.^9+16.*x.^6;

u = fu(xx, yy)';
f = ff(xx, yy)';

%u = reshape(subs(subs(sim_u, x), y), m,n);
%f = reshape(subs(subs(sim_f, x), y), m,n);

g=zeros(m,n);
g(1,:) = u(1,:);
g(:,1) = u(:,1);
g(m,:) = u(m,:);
g(:,n) = u(:,n);

%crop(f)
tic
v = solverPoi(crop(f), g);
toc
surf(abs(crop(u) - v));

abs(crop(u) - v)
